#  This script generates Figures A3, A4, and A5.

rm(list = ls())

library(dplyr)
library(magrittr)
library(ggplot2)


results_judicial <- readRDS(here('regression_results_judicial_const_uncert.rds'))


results_legislative <- readRDS(here('regression_results_legislative_const_uncert.rds'))

results_ind <- readRDS(here('regression_results_indliberties_const_uncert.rds'))

#@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

# Judicial constraints ----

## average all estimates
ave_estimates_judicial <- results_judicial %>% group_by(model_id, term) %>% summarise(average_estimates = mean(estimate), count = n())

# merge back to the results
results_judicial <- left_join(results_judicial, ave_estimates_judicial, by = c('model_id', 'term'))

results_judicial$point_estimate_diff_squared <- (results_judicial$estimate - results_judicial$average_estimates)^2
  
across_vars_judicial <- results_judicial %>% group_by(model_id, term) %>% summarise(sample_across_var = sum(point_estimate_diff_squared) / (n()- 1))


variance_means_judicial <- results_judicial %>% group_by(model_id, term) %>% summarise(mean_variance = mean(std.error^2))

variance_means_judicial <- left_join(variance_means_judicial, across_vars_judicial, by = c('model_id', 'term'))


variance_means_judicial$variance_final <- variance_means_judicial$mean_variance + variance_means_judicial$sample_across_var*(1 + 1/899)

variance_means_judicial$se_final <- sqrt(variance_means_judicial$variance_final)

ave_estimates_judicial <- left_join(ave_estimates_judicial, variance_means_judicial, by = c('model_id', 'term'))

ave_estimates_judicial$lo95 <- ave_estimates_judicial$average_estimates - 1.96*ave_estimates_judicial$se_final

ave_estimates_judicial$hi95 <- ave_estimates_judicial$average_estimates + 1.96*ave_estimates_judicial$se_final

ave_estimates_judicial$model_id <- as.numeric(ave_estimates_judicial$model_id)

judicial_results <- ggplot(subset(ave_estimates_judicial, term == 'v2x_jucon'), aes(x = as.factor(model_id), y = average_estimates, ymin = lo95, ymax = hi95)) + geom_pointrange() + theme_bw() +  geom_hline(yintercept = 0, color = 'red', linetype = 'dashed') + scale_x_discrete(breaks = unique(ave_estimates_judicial$model_id)) + labs(x = "Models", y = "Judicial constraints \n (with measurement uncertainty)")


ggsave(judicial_results, file = here('figure_A3.pdf'), width = 4, height = 4)

# Legislative constraints ----

ave_estimates_legislative <- results_legislative %>% group_by(model_id, term) %>% summarise(average_estimates = mean(estimate), count = n())

# merge back to the results
results_legislative <- left_join(results_legislative, ave_estimates_legislative, by = c('model_id', 'term'))

results_legislative$point_estimate_diff_squared <- (results_legislative$estimate - results_legislative$average_estimates)^2

across_vars_legislative <- results_legislative %>% group_by(model_id, term) %>% summarise(sample_across_var = sum(point_estimate_diff_squared) / (n()- 1))


variance_means_legislative <- results_legislative %>% group_by(model_id, term) %>% summarise(mean_variance = mean(std.error^2))

variance_means_legislative <- left_join(variance_means_legislative, across_vars_legislative, by = c('model_id', 'term'))


variance_means_legislative$variance_final <- variance_means_legislative$mean_variance + variance_means_legislative$sample_across_var*(1 + 1/899)

variance_means_legislative$se_final <- sqrt(variance_means_legislative$variance_final)

ave_estimates_legislative <- left_join(ave_estimates_legislative, variance_means_legislative, by = c('model_id', 'term'))

ave_estimates_legislative$lo95 <- ave_estimates_legislative$average_estimates - 1.96*ave_estimates_legislative$se_final

ave_estimates_legislative$hi95 <- ave_estimates_legislative$average_estimates + 1.96*ave_estimates_legislative$se_final

ave_estimates_legislative$model_id <- as.numeric(ave_estimates_legislative$model_id)

legislative_results <- ggplot(subset(ave_estimates_legislative, term == 'v2xlg_legcon'), aes(x = as.factor(model_id), y = average_estimates, ymin = lo95, ymax = hi95)) + geom_pointrange() + theme_bw() +  geom_hline(yintercept = 0, color = 'red', linetype = 'dashed') + scale_x_discrete(breaks = unique(ave_estimates_legislative$model_id)) + labs(x = "Models", y = "Legislative constraints \n (with measurement uncertainty)")

ggsave(legislative_results, file = here('figure_A4.pdf'), width = 4, height = 4)


# Ind. liberties ---- 

ave_estimates_ind <- results_ind %>% group_by(model_id, term) %>% summarise(average_estimates = mean(estimate), count = n())

# merge back to the results
results_ind <- left_join(results_ind, ave_estimates_ind, by = c('model_id', 'term'))

results_ind$point_estimate_diff_squared <- (results_ind$estimate - results_ind$average_estimates)^2

across_vars_ind <- results_ind %>% group_by(model_id, term) %>% summarise(sample_across_var = sum(point_estimate_diff_squared) / (n()- 1))


variance_means_ind <- results_ind %>% group_by(model_id, term) %>% summarise(mean_variance = mean(std.error^2))

variance_means_ind <- left_join(variance_means_ind, across_vars_ind, by = c('model_id', 'term'))


variance_means_ind$variance_final <- variance_means_ind$mean_variance + variance_means_ind$sample_across_var*(1 + 1/899)

variance_means_ind$se_final <- sqrt(variance_means_ind$variance_final)

ave_estimates_ind <- left_join(ave_estimates_ind, variance_means_ind, by = c('model_id', 'term'))

ave_estimates_ind$lo95 <- ave_estimates_ind$average_estimates - 1.96*ave_estimates_ind$se_final

ave_estimates_ind$hi95 <- ave_estimates_ind$average_estimates + 1.96*ave_estimates_ind$se_final

ave_estimates_ind$model_id <- as.numeric(ave_estimates_ind$model_id)

indlib_results <- ggplot(subset(ave_estimates_ind, term == 'v2xcl_rol'), aes(x = as.factor(model_id), y = average_estimates, ymin = lo95, ymax = hi95)) + geom_pointrange() + theme_bw() +  geom_hline(yintercept = 0, color = 'red', linetype = 'dashed') + scale_x_discrete(breaks = unique(ave_estimates_legislative$model_id)) + labs(x = "Models", y = "Individual liberty \n (with measurement uncertainty)")

ggsave(indlib_results, file = here('figure_A5.pdf'), width = 4, height = 4)


